Center for Turbulence Research 
Annual Research Briefs 1995 


391 


J > 


/ 




7 


Some progress in large-eddy simulation 
using the 3-D vortex particle method 

By G. S. Winckeimans 

1. Summary of motivation, accomplishments, and future plans 

This two-month visit at CTR was devoted to investigating possibilities in LES 
modeling in the context of the 3-D vortex particle method (= vortex element method, 
VEM) for unbounded flows. A dedicated code was developed for that purpose. Al- 
though 0(N 2 ) and thus slow, it offers the advantage that it can easily be modified 
to try out many ideas on problems involving up to N « 10 4 particles. Energy 
spectrums (which require 0(N 2 ) operations per wavenumber) are also computed. 
Progress was realized in the following areas: particle redistribution schemes, relax- 
ation schemes to maintain the solenoidal condition on the particle vorticity field, 
simple LES models and their VEM extension, possible new avenues in LES. Model 
problems that involve strong interaction between vortex tubes were computed, to- 
gether with diagnostics: total vorticity, linear and angular impulse, energy and 
energy spectrum, enstrophy. More work is needed, however, especially regarding 
relaxation schemes and further validation and development of LES models for VEM. 
Finally, what works well will eventually have to be incorporated into the fast parallel 
tree code. 

2 . The 3-D VEM method 

We use the 3-D regularized vortex particle method (=vortex element method, 

VEM) as in Winckeimans & Leonard (1993). The particle representation of the 
vorticity field is then taken as 


£<r(x,t) 



II* - **(*)!! 

a 



a) 


with 7*(f) — u;*(t)vol* the particle strength, £ the regularization function, and a 
the core size. All particles have the same core size, and it remains constant in 
time. Particles usually have the same volume of fluid, vol, associated with them 
(e.g., vol = h 3 for particles initially on an h x h x h lattice). Sometimes however, 
the discretization of an initial condition (such as a torus for discretizing a vortex 
ring) leads to particle volumes that are not quite identical, see e.g., Winckeimans 
& Leonard (1993). Since the flow is incompressible, the particle volume remains 
constant in time. We also define the singular (delta-function) particle representation 
of the vorticity field as 


w(x, t) = 6 ( x - *'(*)) 7*(0 • ■ (2) 

8 
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The velocity field, u*, is computed from the particle representation of the vorticity 
field as the curl of the vector potential, which solves V 2 t/v = — Cb a . Hence it is 
divergence-free. 

Vortex elements are convected by the local velocity 

^x*(t) = u,(x*(<),t) , (3) 

and their strength is subjected to 3-D stretching by the local velocity gradient. The 
general mixed scheme is obtained as (Winckelmans 1989, Winckelmans & Leonard 
1988, 1989, 1993), 

= («Vu,(x'(i),t) + (1 - a) (Vu,(x«(!),!)) T ) ■ T*(0 • (4) 

Three different cases are: a = 1 for the classical scheme, a = 0 for the transpose 
scheme, and a = 0.5 for the symmetric scheme. 

For the present version of the VEM code, Gaussian smoothing is used (Leonard 
1985, Winckelmans 1989, Winckelmans & Leonard 1993): 


/ 1 \ 3/! £. 
(W=y e-V, 

(5a) 

°w-sb-‘(;a)' 

(56) 

K{ P ) = ±{G{p)-ap)) , 
r 

(5c) 

F(p)=±(3KW-(W) , 

(M) 


with C the vorticity smoothing function, G the Green’s function for the vector 

potential (= streamfunction), K the Biot-Savart function for the velocity evaluation, 

F a function used in evaluating the velocity gradient, and p = r/o the dimensionless 

distance. This choice leads to a second order method, provided 0 < h/cr < 1. 

2 

The error function erf(x) is computed using e" x and Eq. 7.1.26 in Abramowitz 
and Stegun (1972). For small p, Taylor series expansions are used to evaluate G, A", 
and F. Notice that, in general, switching from / = f a if £ < #o to / = /& if x > xo 
is programmed without an “if’ statement by making use of a Heaviside function: 

/ = /• + (/»-/«) ^(1 +sign(l,® -x 0 )) . (6) 

With the particle strength exchange scheme for viscous diffusion (Mas-Gallic 
1987, Degond &: Mas-Gallic 1989), we have: 



+ ^ S ( vol? ?'(<) “ vo1 * 7 9 (*)) 
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where »?(/>) = — Note that the Gaussian smoothing is the only one for which 

rj(p) = C(p)- (It is also the natural kernel for the diffusion equation (Winckelmans 
& Leonard 1993).) For non-uniform diffusion coefficients (such as in LES), the 
formulation simply becomes: 



+ ^2 X + V ( x9 ( t ))) ( vol? 7*(<) “ vo1 * 7 ? (<)) 

S 


1 /||X*(<)-X9(<)|| 

— * — 


( 8 ) 


2.1 Particle redistribution schemes 

One needs to maintain the condition that particle cores overlap. In some cases, 
this calls for a particle redistribution scheme. The high order A2 scheme used by 
Koumoutsakos (1993) and Koumoutsakos & Leonard (1992, 1995) was adopted. 
It consists of replacing the whole set of vortex particles by a new set. The new 
particles are located on an h x h x h lattice (hence ail particles have vol = /i 3 ). 
Consider first the normalized 1-D problem with unit spacing. Then, in the A2(x) 
scheme, an old particle located at — | < x < | gives — |x(l — x) of its strength to 
the new particle located at —1, (1 — x)(l + x) of its strength to the new particle 
located at 0, and ^x(l + x) to the new particle located at 1. This scheme is such 
that: 

x" = (-1)" (-|*(1 - x)j + (0)" ((1 - x)(l + x)) + (1)" Qx(l + x)) (9) 

for n = 0,1,2. In 3-D, one applies the scheme as A2(x) A. 2 (y) A2(z). This scheme 
then conserves exactly total vorticity, 1 1 = f v u>dx = 2 , linear impulse, I = 
5 f v x x w dx = 5 Yls x * x 7*) angular impulse, A = j f v x x (x x u>) dx = 
| x s x (x* x 7*). It usually performs very well on energy conservation and well 
on enstrophy conservation. 

Notice that a simpler scheme is the Ai scheme: in that case, an old particle 
located at — | < x < \ gives | — x of its strength to the new particle located at 
— and \ + x to the new particle located at This scheme is such that: 

for n = 0,1. Again, in 3-D, one applies the scheme as Ai(x) Ai(y) Ai(z). This 
scheme then conserves exactly total vorticity and linear impulse. It does not con- 
serve angular impulse. It usually performs poorly on energy conservation and very 
poorly on enstrophy conservation. We do not recommend its use. 

The A2 scheme has been incorporated in the fast 3-D parallel tree code as well 
(Winckelmans et al. 1995). Particle redistribution is programmed using the tree 
code data structure. It runs very efficiently. Its cost is much less than the cost 
associated with the field evaluation. 
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2.2 Relaxation schemes for the particle vorticity field 

The particle representation of the vorticity field, does not constitute a gener- 
ally divergence-free basis (Saffman & Meiron 1986, Winckelmans & Leonard 1988, 
1993). Thus, although the initial particle discretization of a vorticity field can be 
made very near divergence-free, this condition does not necessarily remain satisfied 
in long time computations. A relaxation scheme can be applied, if and when neces- 
sary, which ensures that the particle field, &*, remains a good representation of the 
true divergence-free vorticity field, uv = V x u*. Different approaches have been 
proposed (Winckelmans 1989, Pedrizzetti 1992, Winckelmans & Leonard 1993). 

Notice first that, once computed, the velocity gradient tensor, Vu ff , contains 
all the necessary components to evaluate the true vorticity field at the particle 
locations. This vorticity field is then used in both relaxation schemes considered 
here. Notice also that = Vxu ? = Vx(Vx Vv) = -V 2 Vv + V(V-^ ff ). Recalling 
that V 2 ^ = —£<r, it follows that V(V • Vv) =w, - Gj„. 

The P-relaxation scheme (Pedrizzetti 1992) was developed in the framework of 
singular vortex particles. It is modified to be used in the context of regularized 
vortex particles. At every time step, the particle strength vector is modified using 
the filtering: 

7.».. = d — W+/«^{jM <»» 

where uj a (x g ) is the true local vorticity field and where / is a frequency factor. The 
time scale l/f must be “tuned” with respect to the time scale(s) of the physical 
phenomena under study to give satisfactory results. This relaxation scheme basi- 
cally acts as a “spring” that tries to maintain the particle strength vector aligned 
with the true vorticity vector. This simple scheme is such that: (1) It doesn’t do 
anything to the particle strength vector if that vector is aligned with the vorticity 
vector; (2) It is a simple local operation on the particle strength vector. No system 
of linear equations involving neighbor particles needs to be solved. 

The W-relaxation scheme (Winckelmans 1989, Winckelmans Leonard 1993) is 
based on the functional representation of the vorticity field: one requires that, at 
particle locations, the particle vorticity field be equal to the true vorticity field: 




( 12 ) 


This scheme is best applied after the particle redistribution scheme . The fact that 
the particles are then well-aligned on a regular lattice greatly favors the reconstruc- 
tion of a smooth function from the particle strengths. 

It is also best to use Gaussian smoothing as this smoothing permits a “good- 
quality” reconstruction of a smooth function from the particle strengths in the 
whole range of core overlapping: 0 < h/a < 1.5 . With other smoothings, the 
window of acceptable h/cr is much narrower. 

The W-scheme amounts to solving a system of linear equations involving only 
near neighbors. This is done using an iterative method such as Relaxed- Jacobi (in 
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the tree code) or Relaxed Gauss-Seidel. Notice that the matrix is not diagonally 
dominant. Actually, with Gaussian smoothing and particles on a regular lattice in 

d-dimension, the diagonal dominance is violated as soon as < §: In 

1-D, this means h/cr < 1.25. In 2-D, h/c r < 1.77, and in 3-D, h/cr < 1.99. Thus: (1) 
The higher the dimension, the worse the non-diagonal dominance; (2) The smaller 
the worse the non-diagonal dominance. Since we operate here in 3-D, and at 
h/a = 0.75 — 1.0 or so (to satisfy the core overlapping condition), we definitely do 
not have diagonal dominance. 

At this point, the efficient iterative solution of this system is still a subject of 
active research (A. Leonard, private communication). There appears to be an “op- 
erating window” of h/cr where, although not diagonally dominant, all eigenvalues of 
the matrix are still real and positive. In that case, iterative solvers (with or without 
preconditioning) can be developed. For instance, it is known that the Gauss-Seidel 
iteration converges for any symmetric, positive definite matrix (Golub and Van 
Loan 1983). The matrix here is symmetric. It is also positive-definite as long as all 
eigenvalues remain real and positive. 

2.3 Time integration 

For time integration, the O ((A*) 2 ) Adams-Bashforth scheme (AB2) is used. 
Since this scheme is not self-starting, an O ((A*) 2 ) Runge-Kutta scheme (RK2) 
is used for the first time step (after the initial condition or after each use of the 
particle redistribution scheme). This approach allows one to maintain second order 
accuracy. Numerical experiments have indeed shown that an O ( A£) Euler scheme 
for the start-up step is simply not acceptable. The RK2 scheme is efficiently pro- 
grammed as follows: Euler predictor, Trapezoidal-rule corrector. 


3. Energy, enstrophy, and their spectrum 

A formulation developed by Leonard (1976 unpublished, private communication) 
(see also Leonard 1985, Shariff et ah 1989), is used to compute the energy spec- 
trum. Although developed in the context of vortex filament methods (for which 
the filament vorticity field is, by construction, equal to the true vorticity field), the 
formulation also applies to vortex particle methods as long as the particle vorticity 
field, uv, remains a good representation of the true vorticity field, If this con- 
dition is violated, then the evaluation of the energy and of its spectrum becomes 
very complex, see e.g., Winckelmans (1989), Winckelmans & Leonard (1993), Kiya 
(1993). 

With Gaussian smoothing, the energy spectrum is finally obtained as 


E(k) = e E(k ) with E(k) 


and the total energy as (Winckelmans Sz Leonard 1993) 
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= kjy*’-* dX= \f v H-** <IX - (H) 

The enstrophy spectrum is S(k) = k 2 E(k) and the total enstrophy is (Winckelmans 
& Leonard 1993) 




7 ? ' 7* 


^ w » dx . 

75 75 


(15) 


Notice that the cost associated with evaluating the energy spectrum is 0(N 2 ) for 
each k. 

A special case is the vortex ring of circulation T and radius R (Leonard 1985). 
In that case we obtain for the energy spectrum of the infinitely thin vortex ring: 
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This complements a result presented in Leonard (1985): 

E{k) = (TRf \ £ J\ [kR (1 - n 2 ) 1/2 ) dfi 


(16) 


(17) 


The spectrum, computed using a particle discretization of the vortex ring, is 
presented in Fig. 1 . For small kR, E(k ) = . For large kR, E(k) asymptotes to 

— (Leonard 1985). The fact that E ~ (kR) 2 for small kR is a consequence of 
the non- zero linear impulse associated with the vortex ring, e.g., see Phillips (1956). 
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FIGURE 1. Energy spectrum of a singular vortex ring: ; (kR) 2 /6: ; 

(kR)- 1 / 2: . 

It is found numerically that for a given wavenumber, k , the spacing, /i, between 
the particles used to discretize a ring only needs to satisfy kh < 5 or so in order for 
the discrete sum, Eq. (13), to correctly capture the exact integral, Eq. (16). This 
is very surprising (and not understood at this time) because the integrand varies 
quite a bit from one particle to the next (1 versus roughly sm ^^ ) ). 

For comparison with the single vortex ring, the spectrum of two opposite rings is 
shown in Fig. 2. 

In that case, the linear impulse is zero and one finds that E ~ ( kR ) 4 . Actually, 
with sufficient symmetry, one can even create a system with E ~ ( kR ) s . This 
was obtained by considering six vortex rings on the surface of a cube, see Fig. 3. 
Finally, we find that all vortex loop configurations considered lead to a spectrum 
E(k) ~ fc” 1 for large k and that this appears to remain so when they evolve in time 
using VEM, inviscid or viscous (including LES), see Section 5. 


4. LES and the possible extension to vortex methods 

We consider turbulent flows away from solid boundaries. We also consider the 
general vorticity formulation (Winckelmans 1989, Winckelmans &: Leonard 1993), 
together with an LES formulation which conserves the zero vorticity divergence 
(Mansour et al 1978): 


D 

Dt 


( dui \d u j\ 5 ( ( du dwj \\ . v 

Wi = (“&- +(1 - + aij (&" - &?)) (18) 




log (e/ (rR) 2 ) ^ I log (e/ ( TRj 2 ) 
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FIGURE 3. Energy spectrum of 
of size S/R = 1.25: ; (kR) s 
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for any a. The symmetric case a = 1/2 leads to SijUJj for the 3-D stretching, with 
Sij the rate-of-strain tensor. 

In the basic LES Smagorinsky’s model, the turbulence eddy- viscosity is taken as 

*turb = (C,h) 2 (2S ij S ij ) 1/2 (19) 

with SijSij = S 2 > 0. Typically, C s lies in the range 0.1 — 0.24 (Rogallo & Moin 
1984, Lesieur et al. 1995). Consider the eigenvalues Aj, A 2 , and A 3 of the rate-of- 
strain tensor, with Ai + A 2 + A 3 = Su = V ■ u = 0. The model then produces an 
eddy-viscosity 

t'turb = (C, h ) 2 (2 (A 2 + A 2 + Xl)) 1 ' 2 . (20) 

We certainly agree with Lesieur et al. (1995) that “this simple eddy- viscosity hy- 
pothesis is extremely arbitrary, and substantial progress in LES might be achieved 
by relaxing this assumption”. For the time being, however, a simple extension 
to particle methods of this eddy-viscosity LES model is considered. Since our a- 
regularization of the vortex particle method is basically a Gaussian filter, it appears 
natural to replace the usual Eulerian grid filter h by the particle core size a (Recall 
that h/cr = 0(1)) and to take: 

t'turb = (C,a) 2 VSijSij) 1 ' 2 • ( 21 ) 

Other simple ways of constructing an LES eddy- viscosity have been proposed, 
e.g., the model based on local enstrophy of Mansour et al. (1978): 

iW> = (Cv>0 2 (^)* (22) 

with u>,u>i = u) 2 > 0 and C v « C, (C„ « 0.2 in Mansour et al. 1978). If we recall 
the vector identity, 

S 2 = + V • ( V • (u u)) (23) 

together with the Euler equations, 

^ + V-(uu)«-VP, (24) 

it appears that, to first order, the two models differ by |u> 2 — S 2 « V 2 P. This is 
an interesting result as it could be used to explain the differences in the behavior 
of these two models depending on the pressure’s Laplacian. Indeed, although \w 2 
and 5 2 are both positive-definite, their difference, V 2 P, can have any sign. 

A third model based on the relative rate of change of local enstrophy due to 3-D 
stretching of vortex lines, 


1 D 


(iOiUJi) = 2 


u>j SjjUj 




u>iu>i Dt 


U>iUi 


(25) 
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could also be constructed, e.g., 

k'turb = (C w h f 2^^- (26) 

This model has the property that it “selects” the eigenvalues used to compute the 
eddy- viscosity according to the relative orientation between the vorticity vector, w, 
and the principal axes (eigen vectors) of the rate-of-strain tensor. Indeed, writing 
the components of the vorticity vector in the system of principal axes as , u> 2 , w 3 ), 
this model becomes 


t'turb 




h ) 2 2 


/ AjWj + A2w| + A3CJ3 \ 




u;? + a;? + cj? 


7 


(27) 


Hence a vorticity- weighted average of the eigenvalues is used to produce the eddy- 
viscosity. This model produces a negative eddy- viscosity in regions where enstrophy 
is decreasing (i.e., where vorticity is compressed). Since this is undesirable, one 
should use (version 1) or max(0,a; t 5 I jW_,) (version 2) instead of u>iSiju>j 

(version 0). 

In axisymmetric strain (Ai = A 2 = -A/2 and A 3 = A), the classical LES model 
gives (C, hf %/3 | A|, regardless of the orientation of the vorticity vector. If vorticity 
is aligned with the direction of highest rate-of-strain, the “selective” model (ver- 
sion 1) gives ( C w h ) 2 2|A|. If vorticity is perpendicular to that direction, it gives 
(C w h ) 2 |A|. Since 1 < \/3 < 2, this result also suggests that using C w = C 3 as a 
first “calibration” for the selective model is a fairly good choice. 

In DNS of the Euler equations, the emergence of flat pancake-like structures 
(“potato chips”) that shrink exponentially in time is also observed, e.g., Brachet 
et al. (1992). In that case, two eigenvalues become exponentially large, Aj 
A (— \ — e*/ r ), A 3 « ! + e*/ T ), while the intermediate eigenvalue, A 2 ^ A, 

remains roughly constant. During this self-similar collapse, it is observed that 
the vorticity tends to remain aligned with the eigenvector corresponding to the 
intermediate eigenvalue. Instabilities similar to those leading to streamwise vortices 
in the context of free shear layers are expected to subsequently concentrate the 
vorticity and produce isolated vortex filaments. Modeling such flows with LES, a 
classical model would produce, during the collapse phase, an exponentially large 
eddy-viscosity (hence kill the collapse phase in its early stages by dissipating the 
energy rapidly) while the selective model would produce a fairly constant eddy- 
viscosity (hence dissipate the energy at the end of the collapse phase). Thus, the 
two models would behave quite differently. 

Finally, mixed-schemes that are a judicious combination of the above models 
could also be considered. Whatever the choice, they would have to be validated 
somehow (e.g., using DNS data), including the determination of the “constants”. 

One interesting question is whether one of the simple models above (or a suitable 
mix of them) can produce better results than what is so far obtained with the 
classical Smagorinsky’s model. 
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Note that the vortex method also has potential for the development of a dynamic 
LES model, in the same spirit as in Germano et al. (1991), Ghosal et al. (1992), 
Moin k Jiminez (1993), Ghosal k Moin (1994), Moin et al. (1994), Ghosal et al. 
(1995). For instance, one could compute the velocity fields and derivatives from 
the particle locations and strengths by using Gaussian smoothing at two levels: 
e.g., a filter of width <7 and a filter of width 2<r. This information could then be 
used to “compute” C s in a way similar to what is done so far with dynamic LES 
in grid methods. One must recall, however, that the vortex method with Gaussian 
smoothing is a second-order method. If dynamic LES requires higher order methods 
(as it may . . . , Ghosal, private communication), it might not be feasible in the 
context of VEM. 

5. Fast and slow VEM codes 

A fast parallel oct-tree code, originally developed for three-dimensional N-body 
gravitational problems (Salmon 1990, Salmon k Warren 1994, Warren k Salmon 
1995) has been modified into a fast N- vortex code for vortex flow computations using 
the vortex particle method combined with the particle strength exchange scheme 
for viscous diffusion, with the A 2 particle redistribution scheme, and with both P- 
and W-relaxation schemes (Salmon, Warren k Winckelmans 1994, Winckelmans et 
al. 1995a, b,c,d). 

Gravitation, VEM, etc. are all 0(N 2 ) in complexity: for each of the N elements, 
find the derivatives of the field induced by all N elements. This is the expensive 
part of the computation. The other tasks (particle strength exchange scheme, par- 
ticle redistribution, etc.) are all fairly local operations and are not computationally 
expensive. The use of fast tree codes in 2-D and 3-D reduces the computing cost 
associated with all evaluations from 0(N 2 ) to something much more tractable: 
0(N log N), or 0(N 1+( ) with e « 1, or even O(N), depending on the complex- 
ity of the implementation. The “big-O” notation can, however, be misleading for 
practical values of N and desired level of accuracy. In our implementations of the 
VEM, multipole expansions of order p = 2 are used (i.e., monopole + dipole + 
quadrupole). Particular attention is given to ensuring that the error introduced by 
the use of multipole expansion approximations remains below a desired level for all 
evaluations. A run-time parameter, e to i, determines the maximum allowed error 
bound for any particular multipole evaluation. 

The tree code is written entirely in ANSI C and has been ported to several 
parallel and sequential platforms. Problems with N = O(10 4 — 10® ) and beyond are 
computed on parallel supercomputers. Problems with N — O(10 3 — 10 s ) are also 
computed on the “degenerate” parallel case of single processor workstations. 

For the present two-month “exploratory” work at CTR, it was decided to stick 
with a slow 0(N 2 ) VEM code. (Actually, an all-new VEM code was written for that 
purpose.) Recall that computing an energy spectrum is also an 0(N 2 ) operation 
for each wavenumber k anyway. Although this 0(N 2 ) code sets a limit of N ~ 10 4 
on the number of particles (even on a CRAY C90), it provides for an easy and 
convenient way of experimenting with many ideas: different LES models, different 
particle redistribution schemes, different relaxation schemes, etc. 
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6. Some computational results 

6.1 Twelve rings compact vortex system 

We consider a “compact” vortex system which, by construction, has the following 
desirable properties: zero vorticity (els always in 3-D), zero linear impulse, zero 
angular impulse, and zero helicity, H = \ f y u* x u> dx. Initially, it is formed of 
twelve circular vortex rings, each of circulation T = 1: six rings of radius R — 0.6 (38 
sections per ring, with 9 particles per section (1 in the center with circulation T/2, 
and 8 around the center, at a distance r c = 0.123 and with circulation T/16) laid on 
the surface of an outer cube of size 5 = 1 and with self-induced velocity directing 
them towards the cube’s center, and six rings of radius R = 0.3 (19 sections per 
ring, again with 9 particles per section) laid on the surface of an inner cube of size 
5 = 0.5 and with self-induced velocity directing them away from the cube’s center. 
The total number of particles is N = 3078. The spacing between particles along the 
ring is h 0.10. The two cubes share the same center. The outer cube is directed 
along e x = (1,0,0), e y = (0,0,1) and e z = e x x e 

the inner cube is arbitrarily oriented along e x , e y , 

= (°i |)> e z = e* x e' y , e* = e,/||e z || and e y ■■ 

To ensure core overlapping for a long time, a large value of a = 0.20 is used (hence 
h/a « 0.5). The time step is A t — 0.02. The symmetric stretching scheme is used, 
a = 0.5. The LES model of Eq. (21) is used, with C 9 = 0.1. The W-relaxation 
scheme is used every 10 time steps (with 50 Gauss-Seidel iterations). 

Initially, the energy is E = 1.428 and the enstrophy 5 = 46.21. Following classical 
definition of (isotropic) turbulence, the integral length scale is obtained as 


y . To break the symmetry, 
e z , with e x = (±, 

x e. . 


L = 


J^k-^Ejk)dk 


= 0.490 


(28) 


and the Taylor microscale els 


( £\ 1/2 

= (5-) =0.393 


(29) 


At first, a run without particle redistribution is conducted up to t = 4. Contour 
plots are presented in Fig. 4. The histograms of energy and enstrophy are provided 
in Fig. 5. 

The energy decays due to LES diffusion. Due to vortex stretching, the enstrophy 
first increELses. It then decreases due to vortex reconnection by viscous diffusion. 
Notice that two enstrophy curves are presented. The 5-curve refers to enstrophy as 
defined by Eq. (15). The 56-curve refers to enstrophy defined as 


£b = \ I Cbdx . (30) 

2 J v 

As long els the pELrticle vorticity field, a;*, remains a good representation of the 
divergence-free field, uj (Ti the two curves remain identicsd. Their difference is thus a 
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FIGURE 4. Twelve rings interaction. 3-D contour plots of — 2.0 at t = 0.0 and 
2.0 . 
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FIGURE 5. Twelve rings interaction. : without redistribution: E and £: • , 

£y. x ; : with Ai redistribution: E and £: * (inverted), £y. x . 



FIGURE 6. Twelve rings interaction. Energy spectrums: t = 0: , t = 1: 

, t = 2: , t = 3: , t = 4: ; Jt 8 , it 6 , jfc 4 , it 2 , Jfc" 1 , Ar 5 / 3 : 


global indication of problems with d;* ^ u^. In the present case, it is seen that the 
W-relaxation scheme does a fairly good job at keeping uj a « u;* up to t « 2 or so. 

Energy spectrums axe provided in Fig. 6. It is seen that the high end of the 
spectrum starts filling up at t « 2 or so. This is also indicative of problems with 
This is confirmed by a close look, for all particles, at the amplitude of 
tie and and at their relative orientation. It is also seen that the low end of the 
spectrum does not remain well-behaved as time evolves. The behavior is physically 
acceptable as long as it remains above (kS) 4 . The fact that it evolves to (kS) 2 
indicates that spurious creation of linear impulse has occurred. This is confirmed 
by a close look at the histogram of I(*).. Finally, total vorticity, fl(t ), also does not 
remain zero as it should. This could be somewhat improved by using the transpose 
scheme, a = 0, instead of the symmetric scheme (Winckelmans 1989, Winckelmans 
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FIGURE 7. Twelve rings interaction. Unfiltered energy spectrums: t = 0: , 

t = 1: , t = 2: , t = 3: , < = 4: ; fc 8 , A: 6 , fc 4 , A: 2 , k~ l : 


& Leonard 1993). The W-relaxation scheme, however, does not conserve ft. Again, 
all the above “symptoms” are indicative of problems with Cj ff / 

For comparison, a run with Ai particle redistribution every 50 time steps (and 
with h — 0.10) is also carried out. From N — 3078 at t = 0, this leads to N = 6590 
at t - 1, and to N = 11160 at t = 2. Because of the 0(N 2 ) code, the computation 
cannot be carried out much further than t = 2, and is ended at t = 2.6. Histograms 
of energy and enstrophy axe provided in Fig. 5. The conclusion is that the Aj 
scheme is definitely not acceptable: it dissipates too much energy and enstrophy. 
In particular, it totally overshadows the amount of energy dissipated by the LES 
model. Another interesting result is that the correspondence = ut# is better 
maintained with particle redistribution than without. This confirms that the W- 
scheme is indeed best applied when combined with redistribution. Energy spectrums 
are also provided in Fig. 6. This time, the high end of the spectrum is still fine 
at t = 2. So far, the low end of the spectrum also behaves fine. Although the Aj 
scheme exactly conserves ft and I, it is likely that spurious creation of ft and I 
will also occur eventually due to the W-relaxation scheme and to the symmetric 
stretching scheme. 

One conclusion so far is the following: If one is to do controlled LES with the 
VEM, it must be that the energy dissipation due to redistribution or relaxation is 
less than the one due to LES. A good run might be to use the A 2 scheme every 10 
or 20 steps. (This scheme indeed conserves much better energy and enstrophy, see 
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below.) This could not be done with the present 0(N 2 ) code, however, due to the 
large increase in the number of particles required. 

Another conclusion is that the W-scheme does not conserve SI and I. (Neither 
does the P-scheme.) One further improvement would be to develop a relaxation 
scheme which conserves SI (and, if possible, also conserves I). 

One question arises regarding the “inertial” range of such vortex tubes interac- 
tions. Is there a (A;S)“ 5 / 3 Kolmogorov range that develops? In Kiya (1993), it 
is argued that yes, there is. We claim that no, there is not. In considering the 
filtered energy spectrum, E(kS ), of Fig. 6, it is hard to tell whether there is a Kol- 
mogorov range or not. One finds the answer by considering instead the unfiltered 
energy spectrum, jE(A: 5), of Fig. 7. Then, there is a clear indication that (1) the 
computation blows up (see comments above), and (2) as long as it doesn’t blow 
up, the spectrum remains as ( kS)~ l . This point will become clearer below, on a 
computation that replicates the one presented in Kiya. 

6.2 Six rings compact vortex system 

We consider next another compact vortex system with zero vorticity, zero linear 
impulse, zero angular impulse, and zero helicity. Initially, it is formed of six vortex 
rings, each of circulation T = 1 and of radius R = 0.6 (38 sections per ring, with 9 
particles per section (1 in the center with circulation T/2, and 8 around the center, 
at a distance r c = 0.123 and with circulation r/16) laid on the surface of an outer 
cube of size 5=1 and with self-induced velocity directing them towards the cube’s 
center. The rings are elliptical (in order to break the symmetry) with ab = R 2 
and a/b = 1.25 (top), 0.80 (left), 1.33 (bottom), 0.75 (right), 0.85 (front) and 0.90 
(back). The total number of particles is N = 2052. The spacing between particles 
along the ring is h « 0.10. 

A value of a = 0.14142 « \flh is used (hence h/cr « 0.707). The time step is 
A t = 0.025 and the computations are carried out up to t = 4. The symmetric 
scheme is used, a — 0.5. The LES model of Eq. (21) is used, with C s = 0.2. The 
W-relaxation scheme is used every 10 time steps (with 50 Gauss-Seidel iterations). 

Initially, the energy is E = 1.745 and the enstrophy £ = 65.39 (hence A = 0.365). 
Notice that the application of the Ai scheme to that perfectly fine initial condition 
leads to E — 1.642 (loss of 6%) and £ = 57.91 (loss of 11%) . For comparison the 
application of the A 2 scheme leads to E = 1.741 (loss of 0.24%) and £ = 64.84 (loss 
of 0.83%). This illustrates the superiority of the A 2 scheme over the Ai scheme, 
regardless of the time evolution of the vortex system. 

Three runs were done: one without particle redistribution, one with A 2 redistri- 
bution at t = 2, and one with Ai redistribution at t = 2. Contour plots for the first 
run are presented in Fig. 8. 

The histograms of energy and enstrophy are provided in Fig. 9. Again, the 
energy decays due to LES diffusion. Due to vortex stretching, the enstrophy first 
increases. It then decreases due to vortex reconnection by viscous diffusion. As 
long as remains close to u; a (here up to t « 2), the two curves, £ and £t> remain 
identical. The Ai scheme is again clearly unacceptable. The A 2 scheme performs 
much better. However, it is believed that it should have been used more often (i.e., 
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FIGURE 10. Six rings interaction. Energy spectrums: t = 0: , t = 1: , 

t = 2: , t = 3: , t = 4: ; k 8 , k 6 , k\ k~\ k ~ 5 / 3 : . 


every 10 or 20 steps instead of every 80 steps) to give a better performance. Again, 
this could not be done, due to the 0(N 2 ) computational cost of the code. Finally, 
the correspondence uv = u) a (and hence £ = £&) is again better maintained with 
redistribution than without. 

The energy spectrums are provided in Fig. 10. The high end of the spectrum 
starts filling up at t & 2. This is again indicative of problems with u> a ^ and is 
confirmed by a close look at both £j a and for all particles. The low end of the 
spectrum remains well-behaved as time evolves, with very little spurious creation 
of linear impulse and of total vorticity. The six rings interaction here constitutes a 
“gentler” problem than the previous twelve rings interaction. 

Again, regarding the “inertial” range of these vortex tubes interactions, it is 
again closer to a ( kS)~ l behavior (filtered by the Gaussian) than to a (&S)“ 5 / 3 
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Kolmogorov behavior. 

6.3 Six thin rings inviscid vortex system 

To settle the issue, a run that replicates Kiya (1993) is also done. In that case, 
six circular rings of radius R = 1 and of circulation T = 1 are laid on the surface 
of the cube of size S = 1.25. Each ring is discretized using a single line of 256 
particles (hence h « 0.0245. In Kiya, the high order algebraic smoothing is used, 
with < 7 * = 0.10. Recalling that the self-induced velocity of such a ring is obtained 
as (Leonard 1985, Winckelmans 1989) 



whereas the velocity of the ring with Gaussian smoothing is (Leonard 1985, Winck- 
elmans 1989) 

JU4 

47t R [ \ 


U = 


'?) 


- 1.058 


(32) 


the proper scaling requires that our computation be done with a = 0.05724. Thus, 
these are much thinner rings than before. Hence a wider “inertial” range is expected. 

The computation is carried up to t = 1.5, with At = 0.01. Again, the symmetric 
scheme is used, a — 0.5. This is also a simple VEM computation. Hence, no relax- 
ation scheme, and no redistribution scheme. Finally, this is an inviscid computation. 
Hence, no LES. 

The energy spectrums are provided in Fig. 11. As claimed by Kiya, the filtered 
spectrum, E{kR) suggests a ( kR)~ 5 / 3 behavior. This is purely due to the filter, 
however. Indeed, from examining the unfiltered spectrums, E(kR ), it is clear that 
(1) the behavior remains as ( kR ) _1 for a long time (forever?), and (2) the compu- 
tation eventually blows up (as was the case in Kiya). The histograms of energy and 
enstrophy are provided in Fig. 12. From the difference between the curves £ and 
£&, it appears that the computation blows up at t « 1.2. 

In conclusion, it appears that interactions involving only vortex tubes lead to a 
A:" 1 behavior. It may require the interaction between both vortex tubes and vortex 
sheets to obtain a Kolmogorov-like spectrum. A model involving spiral vortices 
(i.e., rolled-up vortex sheets) is presented in Lundgren (1982). 

6.4 DNS of two rings fusion using the fast parallel tree code 

This work was not done while at CTR. It was done in collaboration with Salmon, 
Warren and Leonard (Winckelmans et al. 1995d). It is also presented here in order 
to illustrate the capabilities of the fast parallel VEM code. We consider a high 
resolution DNS of the fusion of two vortex rings: radius R = 1 , circulation T = 1 , 
Gaussian vorticity distribution with or — 0.10, spacing of the two rings center to 
center S — 2.70, angle of each ring w.r.t. vertical = 20 degs. Each ring is discretized 
with 126 sections and 225 particles per section (i.e., 7 layers, see Winckelmans &; 
Leonard 1993). The inter-particle spacing is then h « 0.05. The computations were 
run with At = 0.05, a = 0.0625, a = 0, v = 0.0025 (i.e., Re = T/u — 400) on 
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Figure 11. 

spectrums: t 
fc -5/3. 


Six thin rings inviscid interaction. Filtered and unfiltered energy 
0: , t = 0.5: , t = 1: , t = 1.5: ; fc 8 , k~\ 



FIGURE 12. Six thin rings inviscid interaction. : without redistribution: E 

and £: • , £b '■ x • 


both 32 nodes of the NAS IBM-SP2 and 64 nodes of the Caltech Intel Paragon. 
Initially, there were 56,700 particles (19 CPU secs per step on SP2-32 and 68 on 
Paragon-64). The A 2 particle redistribution scheme with h = 0.05 was used every 
10 time steps. At the end of the run, there were 218,696 particles (87 CPU secs per 
step on SP2-32 and 236 on Paragon-64). The velocity error was roughly 0.0006 for 
the mean over all elements, and 0.0008 for the max. 

It is seen in Fig. 13 that the diffusion scheme, when combined with the high 
order particle redistribution scheme, correctly captures the fusion process: First, the 
energy and enstrophy losses associated with the A 2 scheme are small enough that 
they cannot be seen in the histograms. (They can only slightly be seen when they 
are differentiated numerically.) Second, the normalized energy decay rate remains 
(almost) equal to the enstrophy, as it should. For comparison, a run without particle 
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redistribution was also done, see histograms in Fig. 13. In that case, the energy 
decay rate is clearly incorrect. Finally, the conservation of linear impulse is also 
much improved by the use of the redistribution scheme. Yet, even with particle 
redistribution, linear impulse starts decreasing at t « 4. It is believed that Cb a is 
then beginning to deviate from u ; (7 . 

At this point, we are also experimenting with the two relaxation schemes when 
used in conjunction with the redistribution scheme. Results obtained so far are 
encouraging, yet too preliminary to be reported. 

7. Conclusions 

The VEM method has gone a long way since its early stages: accurate viscous dif- 
fusion, particle redistribution schemes, relaxation schemes for the particle vorticity 
field, fast and accurate field evaluation on both sequential and parallel platforms. 
This work is still in progress. The time has come to start developing LES models 
suitable to VEM. During this two-month visit at CTR, a dedicated 0(N 2 ) LES- 
VEM code was developed. Although slow, this code could be modified rapidly in 
order to experiment with many different schemes and ideas. Energy spectrums 
could also be computed. Some progress was accomplished in the following areas: 
(1) LES models and how to incorporate them into VEM, (2) energy spectrums and 
how to compute them, (3) particle redistribution schemes, (4) relaxation schemes. 
More work is needed, however, especially regarding (1) relaxation schemes and (2) 
further validation and development of LES models for VEM (which also requires 
that they eventually be incorporated into the fast parallel tree code.) 

It is believed that, when combined with recent developments in vortex techniques 
for wall-bounded flows (Pepin 1990, Koumoutsakos 1993, Koumoutsakos &: Leonard 
1992, 1995, Koumoutsakos et ai 1994), a matured and well- developed methodology 
will permit the simulation of 3-D unsteady problems of engineering interest: flow 
past airfoils including vortex wake, and flow past bluff bodies including vortex wake. 
These body /wake computations will require the merging of the VEM code with a 
Boundary Element Method (BEM) in order to determine, at each time step, the 
vorticity flux required at solid boundaries in order to satisfy no-slip. 
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